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ABSTRACT 


K V ' 

The structural and dynamic design has been done for a 10 KW 
VGAVWT to be set up by the Ministry of Non-conventional Energy 
Resources, Govt, of India and BHEL. The section consists of a 
standard 2 cell structure consisting of a D-shaped thick walled 
spar and thin skin. The blade is provided with stepped internal 
taper to reduce the weight and to bring down the position of the 
mass centre to facilitate centrifugal reef control. The stresses 
are determined using an iterative technique which uses the 
Engineer’s theory of bending and the Bredt-Batho theory for 
torsion of multi-cell sections. The Tsai-Wu failure criterion for 
composite materials is used. 

The Gallerkin finite element method is used to find the 
natural frequencies and the modes of vibration of the blade. The 
periodic aerodynamic forces are decomposed into harmonics to 
determine the dominant excitation frequencies. The Campbell 
diagram is constructed to check for resonance. The effect of 
change in hinge location, shear centre position, thickness and 
taper on the frequencies of vibration is studied. This is done to 
check for resonance due to slight change in the construction of 
the VGVAWT blade. 

A method of modelling high capacity VGVAWT blades as 
composite box-beam spar has been discussed, the variational FEM is 
used to study the response of the structure subjected to various 
loading conditions. It has been shown that by a proper selection 
of ply orientation angle and stacking sequence , the deflections 
in the structure can be reduced. 
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CHAPTER 1 

INTRODUCTION 

Winds result from the fact that the earth’s equatorial 
regions receive more energy than the polar regions. This sets up 
large scale convection currents in the atmosphere. Meteorologists 
estimate that about 1% of the incoming solar radiation is 
converted into wind energy. Since the solar energy received by 
the earth in just ten days has an energy content equal to the 
entire fossil fuel reserves of the world , this means that the 
wind resource is extremely large. One percent of the daily wind 
energy input (i.e 0,01% of the incoming solar energy) is 

equivalent to the present daily energy consumption of the world* 
As the fossil fuels are exhaustible , sufficient power will 
not be produced using these conventional sources of energy 

in the near future. Naturally one has to go in for alternative 
sources of energy such as solar energy, wind energy , tidal 
energy and biomass energy etc. It has long been realized that 
wind is a clean and renewable source of energy from which 

electrical and mechanical power may be generated. 

In India The Ministry Of Non -Conventional Energy 
Resources has identified the wind rich areas which have the 
potential of producing 2000 MW of electric power. The ministry 
has started wind energy programmes by setting up wind turbines 
at various places like Tuticor in (Tamilnadu), Mandvi, Okha 

(Gujrat), Pur i ( Orissa ) etc. 
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1.1 Wind Turbines 

On the basis of the axis of rotation wind turbines can 
be classified into two types namely 

1) Horizontal Axis Wind Turbine { HAWT ) 

2) Vertical Axis Wind Turbine ( VAWT ) 

A HAWT requires a yawing system to align it with the 
direction of the wind. On the other hand, a VAWT is able to extract 
energy from the wind independent of the wind direction and thus 
one can avoid the cost and complexity of a yaw orientation 

sub system. Further advantage of the Vertical axis configuration 
is that it eliminates the cyclically varying Stresses which the 

Horizontal Axis Wind Turbine experience due to their own weight 
. These cyclically varying gravity stress become progressively 

more predominant as the overall size of the turbine increases 
and imposes a theoretical limit on the size of the blade that 

can be constructed. Although VAWT blades experience cyclically 
varying aerodynamic loads, the resulting stresses remain 

substantially constant as the turbine size is increased. A 

major drawback of the vertical axis wind turbines is that they 
do not have a good self - starting characteristic and the 
usual low-solidity < less number of blades ) machines will not 
reliably self start. A reliable self -starting capability can be 
provided by increasing the rotor solidity, but the increased 
rotor weight and cost make this approach unattractive, except 



3 


perhaps for smaller machines. Experience over the last decade has 
shown that VAWTs are best suited to applications at the tens of 
Kilowatts level and higher, where self - start is not important. 
These include the commercially dominant applications, where wind 
turbines are used in parallel with either a grid system or some 
auxiliary power source, such as a diesel engine. 

However, HAWT have a long history and a mature 

technology. Most of the wind turbines in India and other 

countries are of the horizontal axis type. The VAWT is a 

relatively new concept. Hence, data and literature available for 
VAWT are very small in number. 

1*1*1 Vertical Axis Wind Turbine 

There are two types of VAWT 

1) Darrieus rotor 

2) Straight blade type 

A typical VAWT s - of the straight blade type consists 

of a cross arm, blades, tower, electric generator, control and 
foundation. Figure 1.1 shows the general arrangement of a straight 
blade type VAWT. It is seen that each blade is divided into two 
halves and has a airfoil cross section. When wind flows over the 
blade it causes aerodynamic forces ( i.e lift and drag ) to act 
on the blade. It is seen from fig 1.2 that the lift force does 
not act radially . Hence it produces a torque about the axis of 


rotation and causes the blade to rotate. Also it 


is seen that 
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each blade of the VAWT interacts effectively with the on coming 
wind over two arcs, one upstream and one downstream of the 
tower, when it is crossing the wind. For the remaining of each 
revolution, when the blade is moving appr ok i mat e 1 y parallel 

to the wind, there is little aerodynamic interaction. The optimum 
solidity is therefore somewhat higher than that for a 

comparable HAWT. 

1*1*2 : Centrifugal Reef Control Mechanism, 

The power produced by a turbine depends on the wind 
speed , Most of the wind turbines are located in coastal areas, 
since they have a high wind power potential throughout the year. 
The coast lines are often prone to cyclones and hence wind 

turbines may experience severe loading. In addition to this, 
the speed of rotation of the turbine also increases. It is 

necessary to maintain a constant speed of rotation to produce 
electric power with constant frequency. To maintain the speed of 
rotation, the blades are made to incline inward with respect to 
the vertical. This process is called reefing. The wind turbines 
which incorporate this feature are called Variable Geometry 
Vertical Axis Wind Turbines (VGVAWT). 

Presently, reefing is achieved by a hydraulic system 
which senses the over speed of the wind turbine and allows the 
blade to reef through a telescopic cross-arm arrangement. However 
the hydraulic system can fail under extreme wind conditions. 
Dr . K, Ghosh < 1990 ) suggested a novel method of reef control namely 
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centrifugal reef control technique. 

Figure 1.3 shows the general arrangement for achieving 
centrifugal reef control. It consists of a telescopic cross - arm 
which supports two half blades at its either ends. The half blades 
are connected to the end of the cross arm through a base plate 
and a hinge. The end of the the cross - arm is connected through 
a spring to a point on the rotor axis. Each half blade is connected 
through a control strut to the sleeve which houses the cross 
arm. 

When the wind speed increases , it causes a increase in 
the lift force acting on the airfoil. This results in the increase 
of the torque and hence in the speed of rotation .When this 
happens, the cross - arm moves radially outward in the sleeve and 
causes the blade to reef. In the reefed position, there is a 
change in the distribution of the lift along the length of the 
blade because of which the effective torque decreases. This brings 
back the speed of rotation to its original value, 

1*2: Structural and Vibration Design of VGVAWT 

The design of a wind turbine consists of two stages 
namely aerodynamic design and the structural design. The 
aerodynamic design deals with the selection of approximate blade 
profile and other aerodynamic parameters to generate the required 
power. These include determination of the radius of rotation of 
the blade , the length of the turbine blade f the chord length of 
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the airfoil, twist of the section if required etc. From this, 
aerodynamic forces viz. lift and drag can be estimated. 

The structural design commences with the determination of 
various aspects of the blade structure. Normally it is a two or 
three cell structure consisting of a thick spar< D - shaped or 
rectangular shaped ) and a thin outer skin. The spar is the major 
load carrying member. The skin is essentially used to give the 
blade the desired aerodynamic shape although it can carry some 
shear stress. The hollow space adjacent to the spar is filled with 
foam to maintain the shape. 

At present three main materials have emerged as practical 
alternatives for blade construction. These are Wood / Epoxy 
laminates, Glass Reinforced Polyester (GRP), and steel. The use of 
composite materials has the advantage that for the same strength 
the mass of the blade is very less compared to those fabricated 
from conventional materials. Also composite materials make it 
possible to form complex shapes. 

The forces acting on the blade are 

1) The drag and the lift forces (Aerodynamic forces) 

2) The centrifugal force 

3) The gravitational force (weight) 

Since the resultant aerodynamic force does not pass through 


the shear centre , 

it causes both 

bend i ng 

and 

twisting of 

the 

section. For a VAWT 

the centrifugal 

force is 

a 

transverse 

force 

and therefore it 

also causes 

bend i ng 

and 

twisting. 

The 


gravitational force acts in the axial direction for a VAWT . 
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The .thicknesses of the spar and the skin are chosen so that 
the stresses are less than the critical value at failure for the 
chosen material* Since the aerodynamic forces are periodic in 
nature, it is necessary to check that resonance does not occur at 
the operational speed. Both, the natural frequencies and the 
excitation frequencies vary with the rotational speed. Hence an 
interference diagram needs to be constructed to check for the 
resonance , Due to the asymmetry of the cross-section, the 
coupling between the bending and torsion needs to be taken into 
account while finding the natural frequencies. If the spar 
consists of composite material , there can be additional coupling 
effects. 

Apart from fatigue failure and resonance, the third design 
criterion is aeroelastic instability. Some of the important typos 
of instabilities are flap-lag instability, pitch -lag instability, 
classical flutter etc. For blades made up of composite materials, 
instability may not be a critical design criterion as the 
viscoelastic nature of the matrix is expected to damp out any 
ex ter nal excitation. 

1*3 Literature Survey : 

Stress analysis of blades made up of isotropic materials 
or unidirectional fabric is done using the classical theories of 
bending of unsymmetric cross-sections and torsion of multi cell 
sections. These are described in standard texts on Aircraft 

CKb)j. 

structures like fiegson Wrl . 
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.Ramesh Chandra .et.al. C21 have studied the behaviour 
of thin walled composite beam of rectangular cross section* The 
response of the structure at various points was measured using an 
optical system, and the results were correlated with the predicted 
values from a simple beam analysis as well as refined Finite 
element analysis* They found that symmetric ply lay up results in 
bending - twist coupling whereas antisymmetric lay up causes 
extension - twist coupling* 

Alan D* Stemple and S.W.Lee C31 developed a finite 
element formulation to take into account the warping effect of 
composite beams of any arbitrary cross sections. For the 
fibre-reinforced composite beams, results were compared to those 
from a shell element formulation and a good agreement was found 
between these two formulations* 

Inderjit Chopra and E.C.Smith [41 formulated a direct 
analytical model for predicting the effective elastic stiffness 
and corresponding load deformation of composite box - beams. They 
investigated the the importance of three non-classical structural 
phenomena. It was found that torsion related warping, shear 
deformation and two dimensional elasticity play an important role 
in the structural response of composite beams* 

Ramesh Chadra .et.al. C53 presented an analytical cum 
experimental study for the static response of composite rotor 
blades with elastic couplings. They expanded Vlasov theory to 
analyze composite two cell rotor blades. The essence of this 
theory is the reduction of two - dimensional stress and 
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displacement fields ( associates with plate / shell segments of 
the blade ) to one - dimensional stress and displacements 
identified with the beam. The effect of shear deformation and 
constrained warping were not taken into account. The blades were 
tested under unit tip bending and torsional loads and the 
structural response was measured by using a laser-optical system. 
Good correlation between the theory and experiments was observed. 

Most of the literature on vibration of rotary blades pertain 
to horizontal axis wind turbines. William Carnegie £63 applied the 
variational method to derive the equations of motion of rotating 
cantilever bl ad i ng . 

Ramesh Chandra et.al. [7] conducted a theoretical and 
experimental study of the free vibration characteristics of a 
rotating thin walled composite box. The governing 

equations in generalized displacement were derived using the 
Newtonian approach. It was shown that the bending - twist coupling 
influences the vibration of antisymmetric box beams. Extension - 
twist coupling on the other hand, does not influence the flexural 
- torsion vibration of symmetric box beams s i gni f ican t l y . 

C.P.Filipich and M. B, Rosales £83 studied the vibration of a 
uniform beam .£r© ©— bend in g and torsion of - a 
rotating about an axis parallel to the beam. They 
presented an analytical solution for isotropic material. A 
numerical example was solved to bring out the influence of the 
coupling. They compared their results with those obtained without 
considering bending - torsion coupling. They concluded that 
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flexural - torsional coupling is unavoidable. By neglecting it 
some fundamental frequencies are lost and others are in error, 

1*4 Objective of* the present work : 

The present work deals with the structural and dynamic design 
of a variable geometry vertical axis wind turbine (VGVAWT), This 
has been done for an actual turbine to be set up by the Ministry 
Of Non Conventional Energy Resources and BHEL Hyderabad, for 
producing 10 KW power. 

Prior to this work, the entire aerodynamic design of the 
turbine blade has been completed. It is decided to use straight 
{ untwisted ) and uniform ( externally un-tapered ) blade of 
symmetric airfoil section - NACA 0015, The half - blade length is 
chosen as 3,25m and the chord - length as 0,45m. The radius of 
rotation of the blade is decided to be 4.5m, *F 4 io 

connecting strut is fixed at the mass centre of the half - blade. 

The cross section consists of a D - shaped spar and thin skin 
to give the blade the desired aerodynamic shape. The empty space 
within the spar and between the spar and the skin is filled with 
foam to maintain the shape. The blade is given internal stepped 
taper along its length by varying the spar thickness. This 
helps to reduce the mass of the blade and hence the centrifugal 
force acting on the blade. Moreover, it brings down the location 
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of the mass - centre which is desirable for centrifugal reef 
control mechanism* The blade material is E- glass/epoxy as it is 
easily available in India. It was decided to make both the spar 
and the skin from unidirectional fabric (UDF) as it simplifies the 
stress and vibration analysis. Initially a suitable thickness is 
assumed for the spar and the skin. Then Engineer's theory of 
bending and Bredt - Batho theory of torsion for multi - cell 
structure is used to calculate the normal and shear stresses. The 
failure index is evaluated at each critical point of the most 
stressed sections using the Tsai - Wu theory of failure for 
composite materials. If the failure index does not turn out to be 
sufficiently small, -the procedure is repeated with higher values of 
thickness till it is achieved. After finding the spar and the skin 
thickness, the shear centre and the torsion constant ar© 
obtained. 

The blade is analyzed as a rotating beam fixed at one end and 
hinged at the point where it is connected to the strut. The 
coupled equations of motions are first derived from Hamiltonian 
principle. The Galerkin finite element method is then used for 
conversion of these equations into an eigen value problem. The 
eigen value problem is solved numerically to obtain the natural 
frequencies and the mode shapes. The periodic aerodynamic force is 
decomposed into Fourier components to find out which harmonics are 
dominant in the excitation. These harmonics are plotted along with 
the variation of the first few natural frequencies (interference 
diagram) to check for resonance at the operating conditions. 
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.In the end, behavior of composite box shaped spar under 
static loading is studied. It was found that without neglecting 
the effects of warping, the displacement and rotations could be in 
error by as much as 200%, Due to lack of time, stress and 
vibration analysis of composite box shaped beams (with warping 
included) could not be undertaken. 

Several assumptions have been made here for the sake of 
simplicity. The effects of shear deformation and constrained 
warping have been neglected in both the structural and 
vibration analysis. The effect of rotary inertia and the stiffness 
of other structures such as cross - arm, tower etc, have not 
been taken into account in vibration studies. Stability studies 
for the blade have not been carried out, 

1*5 Plan of the thesis 

Chapter 2 deals with structural design of the blade 
considering it to be made of ap uni - directional fabric. The 
engineer's theory of bending and Bredt - Batho theory of torsion 
of multi - cell sections is are used to determine the spar and 
skin thickness, The Tsai - Wu theory of failure for composite 
material is used to determine the failure index. The shear 
centre and the torsion constant of the blade section are also 
aval uated. 

Chapter 3 deals with dynamic aspects of blade design. The 
Galerkin finite element method is used to obtain the natural 
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frequencies and the mode shapes. The aerodynamic force is 

decomposed into a Fourier components to determine the dominant 
excitation frequencies. The first few natural frequencies and 

the dominant excitation frequencies are plotted against the speed 
of rotation to get an interference diagram to check for resonance 
at operational speed. Parametric studies of the natural 
frequencies is also carried out with shear centre location, spar 
thickness, hinge location and taper as parameters. 

In chapter 4 the spar is modelled as a composite box 

beam. Here the shear deformation is included but not the 
constrained warping. The variational finite element method is used 
to study the effects of anisotropic nature of composite materials 
on static response of the spar. This is done for various cases of 

5 

ply layups and loadings. Chapter t gives the summary of results and 

A 

conclusions and suggests scope for future work. 
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Fig l.i ; 


A General Arrangement of 


a Typical VGAVWT 
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azimuthal angle * 


with 
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4 Telescopic cross arm 


5 Damper 


Fig 1.3 : Details of Centrifugal Reef Control Mechanism 
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CHAPTER 2 

STRUCTURAL DESIGN OF VGVAWT BLADE 

This chapter deals with the determination of wall thickness 
and other cross-sectional properties of the blade treating it as a 
beam* The beam is assumed fixed at one end (near the base plate) 
and hinged at the point where it is connected to the strut. The 
stresses in the blade are determined using the engineer’s theory 
of bending ( for unsymmetric sections* ) and Bredt - Batho theory 
of torsion for multi- ceil sections* Since composite material 
( unidirectional fabric ) is used in constructing the blade, 
Tsai - Wu failure theory is used to determine the failure index. 

2* 1 Blade structure 

Based on the aerodynamic considerations, symmetric 
NACA - 0015 airfoil is selected for the blade cross - section. It 
is decided that the blade should be straight ( untwisted ) and 
uniform ( externally untapered ). The chord - length of the 
airfoil and the length of the blade are also selected on the basis 
of power considerations* 

Generally the blade cross - section is a multi - cell 
structure.lt consists of basically two parts the spar and the 
skin. The spar is the main load carrying member. The skin is the 
portion which forms the outer covering and gives the necessary 
aerodynamic shape to the structure. The skin, although normally 
not designed to carry any significant load, can carry some shear 
stresses. The hollow region adjacent to the spar can be filled 
with foam. This helps to to maintain the aerodynamic shape. 
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At present three main materials have emerged as practical 
alternatives for blade construction. These are Wood/Epoxy 
laminate, glass-reinforced polyester < GRP ), and steel* The use 
of GRP has the advantage of an established technology and it is 
possible to form complex shapes. The blades made of GRP have high 
strength, low mass and low stiffness. However the tooling costs 

are high. Wood is used as a material for blade construction for 

low capacity turbines. However there is a lack of information on 

the long term material properties of wood and it can be affected 

by moisture. Steel is used for high capacity wind turbines and has 
advantage of low tooling costs, established technology and 
established long term material properties. The disadvantages of 
using steel include the difficulty in forming complex shapes and 
the high mass of the structure. 

The most common forms of blade construction are shown in 
fig 2*1 . Figure 2.1. a shows a typical section from a wood blade. 
The nose is generally formed as a laminated D-spar, The trailing 
edge is a much lighter glass / epoxy and foam sandwich structure. 
Figure 2,l.b shows a structure made of GRP for low capacity 
turbines in which there is no seperate spar, but ribs are provided 
to provide stiffness. Larger GRP blades are manufactured with a 
filament wound spar as shown in fig 2.1.c. The trailing edge can 
be seperate GRP skin or can be wound as an integral part of the 
blade. Figure 2,l.d shows the use of a box - shaped steel spar 
with GRP skin for high capacity turbines. 

The structure chosen for the the wind turbine blade under 
consideration is shown in Fig 2.2. The spar and the skin are to be 
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made from uni-directional fabric of E-glass / Epoxy with a fibre 
volume of 0.6 and density of 2000 Kg/m 3 , The foam is polyurethane 
and has a. density of 230 kg/m 3 . It is assumed that the foam does 
not carry any loads. 

The blade is provided with a taper. The purpose of the 

taper is to reduce the mass of the structure and to bring the 
position of the mass centre of the section nearer to the base of 
the blade. This is desirable in the centrifugal reef control 
mechanism which requires the strut be connected to the C.G of the 
blade. Uniform taper is not possible because maximum bending 
moment occurs at the intermediate support ( i . e at the hinge 
connected to the strut ). External taper is not desirable because 
it would change the aerodynamic loads along the length of the 
blade and reduce the power output. The blade structure has been 
provided stepped internal taper. It is seen from Fig 2.2 that the 
thickness of the spar varies in stages along the length of the 
blade while that of the skin remains constant. 

We choose the origin of the co-ordinate axis at the mass 
centre of the combined section. The y axis is chosen along the 
chord of the airfoil and the z axis is perpendicular to it in the 
plane of the cross-section. The x axis is along the length of the 
blade < Fig 2.2). 

2.2 Loading on the blade: 

The loading on the blade basically falls into three 
categor ies 

1) Aerodynamic force 

2) Inertial force 
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3) Centrifugal force 

The aerodynamic force acts due to the interaction of the 
ambient with the airfoil. It is normally resolved along two 
directions: along the direction of the relative motion! drag ) and 
perpendicular to it ( lift ). For the purpose of stress analysis 
it is convenient to resolve the aerodynamic f orc^ into components 
along the chord! and perpendicular to it ( F n 1* These are 
given by the expressions 


F 

n 


2 


oc 



F 


F t<*> 

where, 



oc 


U 


V 

oc 


2 V 2 C. C . 

OC t 


2 . 1 


F = Normal force , 
n 

F^ = Tangential force, 
p = Azimuthal angle, 

i/ / 9 

a = Free stream air density = 1.224 Kg/m , 

oc 

U = Relative velocity of the blade * n the 

plane of the blade section, 

V = Free stream wind speed = 10 m/s , 
oc 

C = Normal force coefficient, 
n 

= Tangential force coefficient, 

C = Chord length = 0.45m . 

The variation of [ 3 2 , C and C t e.s functions of 

V n t 

oc 

azimutha) angle <p have been determined from the aerodynamic 

analysis C 9 3. From this the normal and the tangential forces can 
be determined. Figure 2.3 shows the variation of these forces over 



one period. 


The inertial and gravitational forces depend on the mass 
distribution of the structure. Since the structure has stepped 
taper, these forces vary from vary from section to section, 
however on each section they are uniform. These forces are not 
known a-priori as they depend on the thickness. For a typical 
value of t = 6mm the magnitude of these forces are shown in Table 
2 , 1 , 


Tab 1 e 2,1: Values of centrifugal and gravitational forces for 


a typical spar thickness. 


Sect ion 

Centrifugal force (N/m) 

Gravitational force (N/m) 


RO 2 C(pA )+ (pA. ) ] 
c f 

C(pA^> +(pA f )]g 

AB 

4235. 82 

103. 34 

BC 

3593. 68 

86. 95 

CD 

2835.57 

67,25 


In the above table, 

R = Radius of rotation of the blade = 4,5m, 

Cl = Speed of rotation = 65 RPM , 

A = Area of cross - section, 
c 

A^ = Area of the region filled with foam , 

Thus the loads on the blade are 

1) Loads in the x ~ z plane : 

Normal aerodynamic force ( F ) 

n 

Centrifugal force 

21 Loads in the x - y plane 
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Tangential aerodynamic load ( F ^ ) 

3) Gravitational load can be assumed to act in x-y or 
x - 2 pi ane. 

Figure 2. 4. a shows the loading in the the the x - z plane 
while fig 2.4.b shows the loading diagram in the x - y plane. The 


aerodynamic 1 oads 

pass through 

the centre 

of 

pressure 

of 

the 

airfoil while the 

centrifugal 

force acts 

at 

the C.G 

of 

the 


section. These forces cause bending and torsion. The gravitational 
forces cause axial stresses. 

From the table 2. 1 and Fig 2.3 we see that the value of the 
centrifugal force ( which is constant for a particular rotational 
speed ) far exceeds that of the cyclic aerodynamic forces. So w© 
feel that static failure criterion with sufficiently large factor 
of safety can be used instead of fatigue failure criterion. 

S. 3 Engineer's theory of' bending* 

The basic assumption of the Engineer's theory of bending 
is that plane sections of the beam remain plan© after deformation 
produced by the loading. It also assumes that the stress 
component a yy and <y zz are zero. Although a composite material is 
used, the fibres are along the blade length. Therefore the 
anisotropy of the blade material does not affect the derivation. 


The expression for the axial stress a at a point A(y,z) 

x x 

of the cross -section is given by 


ct 

XX 



M 

yy 


/ 


4 


. . . 2.2 
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where, 

P 

A 


M 


ax ial force 
area of cross section 
1 

i - i 2 
yz 


I l 
yy zz 


M 


M I 
y yz 

I 

yy 


M 


M - 
y 


yz 


M I 
2 y z 

I 

zz 


I I 
yy zz 


where, 


M = 

Moment 

about 

the y axis 

y 

M = 

Moment 

about 

the z axis 

z 

I and I = 

Moment 

of inertia of 

y z 

the y 

and z 

axis respec 


I = Product of inertia of the cross - section about 

yz 

the y and z axis. 

In the present case, since the y axis is an axis of 

symmetry, I = 0. So the expression for the normal stress reduces 
yz 

to 


M 


M 


<y 


xx 


y + 


1.3 


zz 


yy 


The above theory is applicable to beams possessing any form 
of cross-sections. However, in the analysis of thin walled 
sections some simplifying assumptions can be made as follows. The 
thickness t of a thin walled section is normal ly small compared 
to other cross- sectional dimensions so that stresses can be 
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assumed to be constant over the thickness. In fact stress can be 
calculated at the mid line of the the section wall and this 

represents the average value over the thickness. Furthermore, in 
the calculation of section properties, we neglect squares and 
higher powers of the thickness and assume the section profile to 
be represented by the mid line of its wall. Thus the expression 
for the moment of inertia I and I can be written as follows 

yy z2 

I = X z 2 't ds 

yy ~ 

i , 2 

1 = I y 1 t ds ... 2. 4 

z z 1,2 

where ^ ^ represents the integration along the mid line contour 
of both cell 1 and cell 2. 

2*4 Bredt - Batho theory of torsion for imilti - cell sections: 


Figure 2.5 shows the blade cross section with the shear 
forces acting on it. It is seen that a part of the shear force 
passes through the centre of pressure and a part passes through 
the mass centre of the section. Since they do not pass through the 
shear centre there is a torque. The shear stresses in the 

section are due to both, the shear forces and the torque. 

The equations for shear stress for a general case containing 

aye • 

m sections j^er developed. A typical multi - cell section is shown 

in Fig 2.6. We choose a different co-ordinate system for each 


ceil. A point on the mid line of 
as the origin for the co-ordinate 
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of the cell j in the clockwise direction. The co-ordinate n_. is 
measured in the outward direction along the normal of the mid 
line. 

The shear stress components associated with the cell j are 

denoted by cr . and a .. These are considered positive when they 
xnj xsj K 

act in the positive n ^ and s^. directions respectively* Besides the 
usual assumptions of the Engineers theory of bending, we further 
assume that 

1) The shear stress component is negligible* 

2) The shear stress component cr . is almost uniform over the 

xsj 

thickness of the wall, 

3) The contribution of P to the stress £^^is negligible* 

Since the shear stress a .is assumed to be almost uniform over 

xsj 

the thickness, it is convenient to introduce the shear flow. For a 
cell j, it is defined by 


To 

wal 1 is 


q 


a xsj ^ 


* . 2,5 


find the shear flow an element on the top / bottom 

considered* It implies 

, j 


3 


t xx 


3 


= 0 * 


2 * 6 


Integrating we get 


q 


j 



where 





0 


s . 

J t d<y 

XX 

<f?x 

at the origin of 



. . . 2.7 


co-ord inate s . , 

J 


is the shear flow 
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Differentiating 2,3 and under the assumption 3 above we have 


da 


xx 


V y 
y 


V Z 
z 


#x 


zz 


yy 


2. 8 


Substituting 2,3.d in 2. 3. c gives 


q 


o 


- C 


V Q ' 

__y y 

I 

zz 


V Q 

+ 3 


i 


yy 


2, 9 


where 


s . 

i j 

Q = / y t ds 


s . 

j 

Q 3 = / z t ds 

z 

0 


2 . 10 


As seen from Fig 2.7, at the junction of the vertical wall with 
the top or bottom wall, conservation of shear flow implies the 
following expression for the shear flow in the vertical wall 
common to the cells j and j + 1 


q 


i 


, 3 * 1 


in the s. direction, 
3 


. . . 2,11 


The expression 2,9 and 2.10 for the shear flow involve m 
unknowns q ^ (j = 1, 2, . ,.m ). To find these, we need m 

equations. We get one equation by equating the moment of the shear 
flows q^ about the mass centre G with the torque M^. The 

equation has the form 
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E * 

= I 


^ q ^ 


. . . 2.12 


where is the perpendicular distance of the cell element As . 

I '* 

from G. P * i s considered positive if it is in the positive n. 

' g J 

d i rec t ion. 

To get m - 1 more equations, the deformation of the cross 

section is considered. It is assumed that the cross - section is 
free to deform in the axial direction ( i.e it is free to warp >, 
but within their own plane they can only undergo rigid 
displacement. The most general rigid displacement consists of 
translation along y and z directions and rotation about the x 
axis. However there exists a point in the plane of the cross 
section which remains stationary. This is called the centre of 


twist and it coincides with the shear centre >4 




The displacement components can be expressed in 

terms of the twist about this point. Let V . and V . be the 

sj n j 

tangential and normal displacements of point a point p. We can 
express these as ( see Fig 2,7): 


1) In thetop and bottom walls 


r <p 
s 


r <p 
n 


2 . 1 3 % a 


2) In thevertical wall common to the cells j-1 and j 


V . - V . < = r <p 

s,3 s, 3 - 1 s 


n, j- 1 


2. 13. b 
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3) In the v© r t i ca l wall common to cells j and j + 1 


V . - V 

s > 3 s, j+ 1 


r <b 
s 


V , = - r <p 

n, j + t n r 


2, 13. c 


The shear strain component £ xn ji s assumed negligible. The 

other shear strain component s . is given by 

x s j 


<?u . 

J 

<?s . 
J 


. . . 2.14 


The strain £ . and the shear flow q are related by 

X S 3 


. . . 2. 15 


Using the fact that # du = 0 for j th cell, we get 


-S-r ds . - / -V- ds . - / — -3- 

J 3-1, j 1 * 3 3,3*1 


= 20 . 


j dx 


. . . 2* 16 


where <ti represents integration over the entire ceil j, / 

j - 1 , 3 

and / denote the integration over the region common to j-l and 

j > j + 1 

j th cell and the region common to the j and j+1 cell respectively 


O . ^ 

J 


1 I r s ds 3 

S" I p g ds j 


2, 17 


is the area enclosed by the cell j. 
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We get similar equations for the other cel Is. We note that for the 
first cell, the second integral does not exist, while for the last 
cell, the third integral does not exist. 

How we have ( m + 1) equations 

1) m equations of the type 2.16 
2 > equat ion 2.12 

These are solved to get ( m + 1) unknowns i.e q ^ ( j = 

o 

and j— — . Once q ^ is known we use equation 2.9 to find the 
dx o ^ 

shear flow. For our case m = 2, Therefore there are three unknowns 

, 1 2 , <i<p 

o dx 

2*5 Determination of torsional constant of the cross - section: 


The torsional constant is a property which is essential for 
the vibration analysis of the turbine blade. To determine it, 
consider pure torsion of a multi - cell section as shown in Fig 
2.8, Here since the normal stress zero everywhere, equation 
2.6 implies that 



J 


0 


. , . 2. 18 


This indicates that the shear flow associated with each cell 
is constant. Then equations 2.12 and 2.20 reduce to 


m 

tt t = 2 E q J Cl. 
j = 1 


2. 19 
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ds . 
J 


t 


j- 1 


/ 

j~ 1 > j 


ds . 
3 


t 


q 


j + 1 


/ 


3 , j + 1 



2 ^°j 


d<£ 

dx 


, , , 2,20 


The equations 2.19 and 2.20 comprise ( m+ 1 ) equations in 
{ m+ 1 ) unknowns namely and < j = 1, 2, ... m ). First we 

find q~* in terms of ~~~ by solving equation 2.20 . The 

solution can be expressed in the form 




dx 


2.21 


where oi^are the constants which depend on the shapes and 
thickness of the cells. Substituting 2.21 in 2.19, we get 


M 


t 


jLi J 


dx 


where 


J 


2 E 

j=l 


a . ft . 
J J 


. . . 2. 22 


. . . 2. 23 


This is the expression for torsion constant of the thin walled 
closed multi - cell section. 


2 . 0 Determination of" shear centre: 


Determination of shear centre is essential for 
analysis of the turbine blade. It also provides an 
method to determine the shear stress in the structure. 

Earlier, in the first approach given in section 


v ibrat ion 
a l ter nat e 

2. 4, the 



moment equation was used as one of the equation to determine the 

j 


cons tan ts 


q ~ - Here we use the fact that when the shear forces 


V and act at the shear centre, there is no twisting of the 

cross - section ( see Fig 2.9 ). With this condition, equation 

2. 16 reduces to 


q 3 ds 

2 L 

IJt 


S 

j-1, j 


qj ’ 1 


v t 


for 


f 

j + 1, j 

1 , 2 , . 


q j+1 ds 

3 


(Jt 


m. 


= 0 

2, 24 


For convenience, the expression 2.7 for the shear flow q J can 
be decomposed as 


where, 


q 3 = q 3 (V ) + q 3 (V ) 

y 2 . 


q 3 (V ) = q 3 

y oy 


V Q J 

y y 

i 

zz 


... 2.25 


... 2 . 26 . a 


q j (V ) = q i 

z oz 


V Q ' 
z z 


yy 


2. 26. b 


Here q~* (V ) is the shear flow due to shear force V and q'* <V ) 
y y z 

is due to V . Substituting expression 2. 26. a and 2.26. b 

separately into 2.24, we get 2m equations in 2m unknowns in 2m 

unknowns q ^ and q *"* for (j = 1,2, m ). By solving these 

n oy oz 

equations, we find the shear flows due to V and V from the 

y z 

expression 2. 26. a and 2.26.b. 

To find the position of the shear centre, we equate the the 

moments of the shear flows q^ (V ) and q^(V 1 about the mass 

y z 

centre those of V^and V^, This gives the following expressions for 


the co-ordinates of the shear centre 
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m 


T f q J (V ) P p ds, 

i = if — „ z g j 


(V„) p„ d Sj 


y G — j 


V 

■ y 


. . . 2.27 


Note that for the section under consideration, should come out 
to be zero as y - axis is the axis of symmetry. 

Once the shear centre position is determined, the shear flow 
due to torsion can be determined as follows. Note that the torque 
about the shear centre is given by 


M 


tE 



V 

zz 


...2.28 


where y is the co-ordinate of the centre 
P 

Then using the theory of pure torsion of 

described in the previous section, we ge 

M by solving equations 2.19 and 2.20 
t E 

V and V are given by equations 2, 26, a 

y z 

these shear flows, we get the net shear 


of pressure (Fig 2.5). 
multi - cell sections 
t the shear flow due to 
. The shear flow due to 
and 2.26.b. By adding 
f 1 ow due to V , V and 

y z 


In this work , the shear stresses were found by two methods: 
by applying the method described in section 2.4 .and the method 
mentioned in this section. The values obtained by these methods 
match very well. This shows that the position of the shear centre 


has been determined accurately. 
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2.7 Tsai ~ Wu failure theory 

The aerodynamic forces acting on the blade are periodic in 
nature. Hence fatigue failure criterion should be used for 
determination of the wall thickness. However, reliable values for 
the fatigue strength of the blade material are not available. 
Moreover, the periodic forces are small compared to the 
centrifugal force which is a static force. Hence the static 
failure criterion has been used. 

DoJ 

We use the Tsai-Uu failure theory* in this theory, for plane 

A 

stress problems, the failure index at a point is given by 


FI - F t ♦ F 2 ^ * F u c* . F 22 ov, 2 * F 33 <r 12 2 + 2 F^ ^ ^ 


where, 


... 2. 29 


t 

1 


c 

1 


11 


X A 
t c 


22 


Y . Y 
t c 
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12 


( F ) J F F 
12 * 11 12 


2. 30 


Here, X = Ultimate tensile strength in the longitudinal 

(fibre) direction or direction 1. 

X ^ = Ultimate compressive strength in the fibre direction. 
Y^ = Ultimate tensile strength in the transverse 
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direction or direction 2. 

S = In plane ultimate shear strength, 

o 1 ^ y a 12 = s ^ a ^ e ° ^ stress with respect to directions 

1 and 2. 


The value of the coefficient F,_ is determined from a biaxial 

12 

test. Performing such a. test on a ply is not straightforward, Tsai 

rvS'J 

and Hahn have shown that the value of the coefficient F*_ lies 
> 12 

between the limits of - 1 and 1. In the absence of a measured 

# 

value, F ^ generally set to -0,5 without loss of accuracy. 

According to this theory, failure occurs when the failure 
index exceeds 1. To account for uncertainties in the determination 
of the loads, stresses and strengths, a factor of safety n is 
used. Then failure is said to occur when 


F I = — ... 2.31 

n 

For static case, n can be chosen as 3. However the aerodynamic 
forces a te periodic in nature,. Therefore to avoid fatigue failure 
due to these forces, the factor of safety should be increased. 
Normally it is taken as three times the static value. Thus, the 
value of n = 10 is reasonable, then for failure , FI should be 
1 ess than 0.1. 

Since the blade is made of unidirectional fabric with fibres 

-srt^e.m'&cL ^ 

along the length direction of the blade, the axes 1 and 2 of 
equation 2.29 coincide with the x and s axes. 
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£• 8 Del er mi riat i on of spar /skin thickness 

The spar thickness in section AB i.e t^ ( Fig 2.2 ) is 

determined iteratively as follows. In the first iteration, some 

value of t ^ is chosen. This enables us to calculate the 

centrifugal and gravitational forces. The reactions at the fixed 

and hinged support are then determined using the procedure for 

statically indeterminate problem. Next, the axial force, shear 

forces and bending moments ( P, V , V , M , M ) are obtained at 

y z y z 

the three critical sections E, B and C (Fig 2.4). Then the 

normal stress a at the critical points of the three critical 

xx 

sections is found by applying the theory of section 2.3. Then 

1 2 

the shear flows q and q and a « and o* _ are determined at the 

critical points of the critical sections using the theory 

described in section 2.4. Finally the failure index is 

calculated by using the Tsai - Wu theory of section 2.7. If the 

failure index exceeds 0.1, the next iteration is carried out by 

increasing the thickness. Iterations are continued till the 

failure index becomes less than or equal to 0.1. Once the 

spar/skin thicknesses are determined, the torsion constant J and 
the co-ordinate of the shear centre y^ are calculated using the 
procedures of sections 2.5 and 2.6 . 

A software was developed to perform the structural design. 

The step wise procedure used for this purpose is as follows 
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8. 8 Results: 

The material properties of the material used for the VGVAWT 
blade are given below in Table 2.2 C 10 ]. 


Tab 1 e 2.2: Properties of VGVAWT blade material 


Material used 

E glass/ ©poxy unidirectional 

fibre composite. 

ELASTIC CONSTANTS I Gpa ) 

Modulus of elasticity in 

fibre' direction ( E^) 

45 

Modulus of elasticity in 

the matrix direction ( 1 

12 

Modulus of Rigidity ( fj ) 

5. 5 

STRENGTH PROPERTIES ( Mpa ) 

Ultimate tensile strength in 

fibre direction < X ) 

1020 

Ultimate compressive strength 

in the fibre direction ( X ) 

c 

620 

Ultimate tensile strength in 

the transverse direction (Y^J 

40 

Ultimate compressive strength 

in the transverse direction(Y) 

c 

140 


In plane shear strength (S) 


70 


















Table 2.3 shows the stresses and the failure index ( FI ) 


a 


the most critical point of each of the three critical sections for 

various values of t, , 

1 


Tab l e 2.3: Stresses and failure index at the critical points of 

the VGVAWT blade. 

The results are given for the following ratio of skin/spar 
thickness 


t ^ { mm ) 


Max i mum 
stress 
( Mpa ) 
and 
FI 


<y 



Point E 
spar thick, 
t ( mm ) 


-72, 989 


15.896 


0. 0754 


Point B Point c 

spar thick, spar thick. 


t^(mm) 


- 87,333 


16.737 


0. 0881 


= t { mm ) 
3 


-23. 001 


1 1 . 855 


0. 03439 



a 

XX 

-65. 195 

- 77.972 

-20. 390 

5 

a 

x s 

12, 960 

15.216 

10.004 


FI 

0. 0594 

0.0751 

0. 03437 


<y 

XX 

-60, 426 

-72. 192 

- 18, 662 

6 

Cf 

X s 

11. 104 

14.086 

7.68 


F I 

0.0566 

i 

0.0665 

0.0183 


At the points B and C there is a sudden change in the cross - 
section. This gives rise to a stress concentration at these 
points. There is no standard data available for estimating the 
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stress concentration effects for the present section* We assume 
the stress concentration factor to be 1*5, Now considering the 
effects of stress concentration, failure is said to occur when 

FI < ^4- i.e FI < 0.06666 ... 2.32 
1*5 

we see from the above table that this condition is satisfied at 
the most critical point of section B when the spar thickness 
t ^ = 6mm* For any other value of spar thickness (t^> less than 6mm 
relation 2*32 the FI is grater than that required to prevent 
failure* Hence we choose the spar thickness in section AB (t^ ) as 


6 mm. 




Force ( N/ 
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CHAPTER 3 

VIBRATION ANALYSIS OF V6VAWT BLADE 

This chapter deals with determination of the vibration 
charec ter is t ics of the VGVAWT blade. The equations of motion 
governing the vibrations are developed using the Hamiltonian 
principle. The turbine blade is treated as a beam rotating about an 
axis, which is parallel to the length of the beam. The effects of 
shear deformat ion, rotary inertia and constrained war ping have not 
been taken into account. 

The finite element method is used to convert the equations of 
motion into a generalised matrix eigen value problem. The beam is 
discretized into a number of linear elements, with each node having 
five degrees of freedom. Natural frequencies and modes of vibration 
are determined by solving the eigen value problem. Parametric study 
is carried out to examine the effects of location of shear centre, 
location of hinge support, wall thickness and taper. 

The periodic aerodynamic force acting on the blade is 
decomposed into Fourier components to determine the most dominant 
excitation frequencies. The variation of first few natural 
frequencies as well as the dominant excitation frequencies with 
the rotation speed is plotted as a Campbell diagram to check for 
resonance . 

3*1 Equations of motion for* rotating beam 

The equations of motion for a rotating beam of arbitrary 
cross - section are derived in reference C l.For the sake of 
completeness, this derivation has been reproduced below with 
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certain modifications. Figure 3, 1 shows the cross - section of a 

beam in y-z plane where ( x, y, z ) is an inertial co-ordinate 

system with origin G at the trace of the rotation axis. Normally 

displacements are expressed with respect to a body - fixed 

co-ordinate system which rotates with the beam. Let ( x , y z > be 

1 X > 1 

such a rotating system with origin at the centroid G at the beam. 

In this system x - axis is parallel to x - axis but not 

coincident. Further, ( y , z ) are the principal inertial axes of 

the cross - section. Therefore 1 =0. Note that the co-ordinate 

y l Z i 

system used in chapter 2 is identical to < x^> y , z ) . The 

angular velocity of the system is denoted by ft . 

Let the displacements of the particle which occupies the 

position P ( x^, z^ ) be ( u, v, w) with respect to the 

co-ordinate system ( x^,, y^, z , Then the co-ordinates of the 

displaced position P* become ( x^+ u, y^+ v, w^+ w) . For the sake 

of convenience the shear centre E is assumed in the third quadrant 

of the c-ordinate system ( y , z ). Let { y , z )be the distances 

XL L fc, 

between E and the centroid G parallel to the y ^and z ^respect i ve i y. 

To facilitate the derivation we introduce an auxiliary rotating 

co-ordinate system ( x _ , y z_). Like ( x . , y z 4 >, this also 

Ill 

rotates with ft but it has origin at 0. The axis is coincident 

with x and the y^ axis is so chosen that ft is always along the 

time 

line 0E. Thus y - axis makes a constant { independent Jangle a with 
2 A 

the y .axis. If the y .axis itself can be chosen parallel to the 

* 1 X 

l i,ne 0E, then there is no need to introduce the auxiliary 
co-ordinate system. For a stationary beam under torsion, rotation 
of the beam cross -section takes place about the shear centre E 
(and not about G). Hence the eccentricity of the rotating beam with 
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respect to the axis of rotation is represented by the distance OE 
rather than OG. In our case, QG is the radius of the turbine which 
has been determined from aerodynamic considerations. The distance 
EG has been determined in chapter 2, From this the distance OE can 
be calculated. 

From fig 3.1 the position vector r ^of the point P*can be 
expressed as the sum of four vectors OE % EG GP > and PP# ? In 
the component form, it can be expressed as 

+ r k ... 3, 1 
z 


r 

y 

= V cos { 

ot 

+ at ) 

- w 

sin 

( a. + Ot ) 

+ e cos 

< 

ot 


r 

z 

= V sin ( 

a 

+ at > 

+ w 

cos 

( a + Ot ) 

+ e s in 

< 

ot 


V 

= y E + y l 

+ 

v, 








W 

= Z E + Z 1 


w. 




* 

* - 

3. 

2 


The velocity of the particle P is given by 

r k ...3,3 

z 

where the superscript 'indicates the partial derivative with 
respect to time. Now, the kinetic energy of the beam is given by 

K = — ~ — p S | vel > | ^ dA dx . . . 3.4 

2 voi 

where p is the density of the beam material. Differentiating the 
expression 3.1 with time and substituting it in equation 3.4, 



u 


r j 

y 




we get. 
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L ^2^2^2 2 2 22 

K = / f — p { u + v +w + ft[V+W+e 

0 A A 

+ 2s { V cos a - W sin a. ) ] 

' * « 4 

+ 2Q C -V W+ W V+ e < V sin a + W cos ot >]} dA dx 

« « 4 3 « 5 

where A is the area of cross - section and L is the length of the 
beam* 

Before writing the expression for potential or strain energy, 
we make some simplifying assumptions about the deformation pattern. 
These are the same as those made in chapter 2. Thus we use 
engineer’s theory of bending and Bredt - Batho theory of torsion to 
get the expressions for bending and torsion deformation and 
stresses. Note that while doing so, we are neglecting the effect of 
shear forces on bending deformation. Further we assume that the 
effect of constrained warping on axial stress is small . Then, the 
expressions for the displacements, after neglecting the axial 
deformat ion, become 

u = -v’y-j - w’ Zj, 

7 = V ( K, t I - 6 ( K, t I ( Zg + Z 1 f , 

w = w(x, t ) + 0 ( x, t ) ( I. 

... 3.6 

Here, v and w are the bending deflections i.e. the translations of 
the cross - section along and z^ axis respectively and Q is the 
rotation of the cross - section about the shear centre. Further the 
superscript * represents the partial derivative with respect to 
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x , The expression for the non-zero stresses are given by 


E { v’ ’ Yl + w’ * z ), 


^ * <x t yt e { j = 1, 2 ) 

KS j 


, , 3,7 


where E and are the elastic properties of the blade material, t 

is the wall thickness and c* . , which are constants, are the 

3 

solution of the equation 2,19 and 2,20 of chapter 2, The 
superscript j refers to the ceil number and s is the co-ordinate 
along the mid line of the wall. 

Now, the expression for the strain energy of the beam can be 
written as 


2 2 

U s — / { / C ■ ^ ^xx + — ~ ^xsj 3 dA } dx 

^ Zb zn 

0 A 

, , , 3,8 

Carrying out the integration over A, this expression reduces to 


u = —1— / C E l {v 1 * I 2 + E l (w 11 ) 2 + fj J (0) z ] dx 

2 ^ Vl 

« • , 3,9 

Where, I and l are the moments of inertia about y and z axes 


y i y i 


z i z 1 


‘espectively and J is the torsion constant. 
By Hamilton’s principle, we get 


6/ ( U- K ) dt =* 0 

t. 


, , . 3, 10 


Where, 


,£ v =<&w 


JJy* s £, * 3T jE* t& — 0 at t 3 tp t 2 , , , 3,11 
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After integration by parts with t (once) and with x 
(twice) and using the conditions at t ^and t ? as well as the we get 
the following equations for the bending ~ torsional vibrations of 
the beam 

EI ( v' * • ' ) + pA ( v* ) - pi ( v* ” ) - pA Cl 2 ( v ) 

11 Z 1 Z 1 

-2pA O ( w ) - pA z £ ( & ) - 2pA Oy £ ( 9 ) + pA Q 2 z £ ( & ) = pA0 2 f, 

El ( w” ” ) + pA { w* ) - pi ( w' ’ ’ ! - pA Q 2 ( w ) 

y i y i y i y l 

* • * O O 

+ 2pA a ( v ) + pA y £ ( e ) - 2pA Oz £ ( <9 ) - pA O' y £ ( (9 ) = pAO g , 

pJ ( # ) - pi ( & ) + pi ( Q ) + pA ( z_ v 11 ~ y_ w* ) 

p p E E 

2 2 
- pA O ( z v - y w ) - 2pA ft ( z w + y c v ) = pA ft eh * 

be b b 

. . , 3, 12 

In the above equation 

f = y E + e cos ol , g = z £ - e sin c* , h = z £ cos a. + y g sin c* 
and 

) 

shear centre. 


c V* • * - (p/E ) v* ' ] 

lo ■ o. 

l L 

c v” 6V’] | 0 

= 0 


I k 

1 ^ 


C w’ ’ ’ - (p/E ) if’ ’ 1 

6w | 0 = 0, 

C w* ’ <5w’ 3 | Q 

= 0 

« < 




C 3 <56 

o 

u 

O 

r- i 

• * 

. 3. 13 


I 


I 


I 


A ( y T 


+ z. 


P y l y 1 z t z l E E 

is the polar moment of inertia referred to the 

The boundary conditions obtained are as follows: 
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For the case of the turbine blade under cons Iderat ion, y ^ - axis is 
an axis of symmetry. Hence z = 0. Also due to its small influence, 
rotary inertia terms can be neglected. In other words, terms 
containing the cross derivatives of fourth order can be discarded. 
When this is done, we get the following equations 


El (v MM ) + pA ( v') - pk Cl (v) - 2 pk O { w ) 

z i 2 i 

- 2 pA Oz E ( 6 ) = pA ft 2 f, 


E l 


Y l y l 


( w' ’ ' ' ) + pA ( w* > - pA ft ( w > + 2 pA ft ( v ) 

' ' 9 2 

+ pA y £ { 0 - ft 6 ) = pA ft g. 


pJ ( 0’ ’ ) 


- pi ( 0 ) + pi ft < 0 > 

P P 

2 2 

+ pA y <-w"+2ftv + ft w)= pA ft eh ... 3. 14 

E 


The boundary conditions for these differential equations are 
obtained by deleting the second term from the first and the third 
condition of 3.13 as these are the rotary inertia terms. 


3.2 Formulation of the eigen value problem 


We can use the method of separation of variables to solve the 
system of equations 3.14. Let the solution be of the form 


CENTRAL 

l i T 

<«*#«. A. 


liBRARK 

< AN PUR 

mm 
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V 

= P ( X ) p< t ) 

+ V ( X ) 

o 


w 

= Q < x ) q ( t ) 

+ w ( x ) 
o 


e 

= R ( k ) r ( t ) 

+ & ( X ) 

n 

. . . 3. 15 


where the functions v , w and & represent the solution of the 

o o o 

time - independent inhomogeneous part. The boundary conditions for 

v , w and 6 can be obtained from equations 3.13. 
o o o ^ 

The equations governing the homogeneous part of the solution 
( P < x ) p( t ) ...etc ) can be obtained by substituting 3.15 into 3.14 
and by dividing by p(t), q(t) and r(t). We get the following 
system of equations, 


E I P 
Z 1 Z 1 


* £*££' - p a o 2 p - - 


2pAOy Rr 


= 0 , 


El d”” + 


y l y l 


pAflq. 


2_ + 2pAOP p + pAy E Rr _ pAn Y E Rr 


= 0 , 


pJR 


. . Plj Rr 


-£ + pa 2 1 r 

r p 


pAy dq' pAQ y Qq 2pAfty Pp 

— E + “ E — = o. 


... 3. 16 

Since the free motion should be harmonic, p(t),q(t) and r(t) 
should satisfy the following relations. 


P % = 


~<P 2 P , 


k' 


= “/? 


2 

r = -77 r 


3. 17 


where (3 and 77 are real constants. In order to satisfy the 



equations 3.16, the following relations must hold 



We define the following terms 


P A 


El 


2 1 Z 1 


2 


c 


P A 
/Jj 


. . . 3. 18 


pA 


E l 


y l y l 



> 


X 


2 





3, 13 


By using 3.18 and 3.19 the system of equations 3.16 can be 
written as. 


# ' : , ' 

< < < > 09 9 9 2 

P - a Z ( X + Q )P - 2a DQ - 2a QygR = 0 , 

ii i« o o 9 99 22 2 

Q - b Z ( X Z + CT)Q. - 2b Z X z OP - b ( X^ + Cl ) y £ R = 0. 

r' * + d 2 ( X 2 + o 2 ) R + c 2 y E < x 2 + 0 2 )Q + 2c 2 y E X 2 OP = 0 

... 3. 20 


These equations show that, when the cross - section has both 
the axes of symmetry, the coupl ing is only between the two bending 
modes { flap-wise and edge-wise ) and it is caused by the 
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rotational speed ft* On the other hand, when ft = 0 and the cross 
section is completely unsymmetr ic , there is coupling between two 
bending modes as well as between each ^bending mode and torsion* The 
coupling between two J&$raT bending modes can be eliminated by 
choosing principal axes of inertia ( equation 3*12 )* Additionally 

if there is one axis of symmetry < say y^ ), then coupling between 
the bending about the perpendicular axis ( i.e z ) and torsion 
also can be eliminated ( equation 3,20 ), When both ft and asymmetry 
are present, there is one additional coupling term involving each 
bending mode and torsion. The equations 3.20 can be written 

in the matrix form as 

C A 3( x'”’) + c B It X * ’ ) + C C H X ) - X. 2 c D H X > = {0} 

. , . 3.21 


where, 

[ X T J = C P Q R 3 




2.2 

-aft 

-2a 2 ft 

-2» 2 «y E 

C C 3 = 

0 

-b 2 0 2 

,2.2 
-b ft y £ 


0 

o 2 « 2 y E 

d 2 ft 2 
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C D ] 


2 

a 

2b 2 Q 

-2c 2 y E 0 


o 


b 


2 


-c y t 



. 3 . 22 


The boundary conditions for P, Q and R can be obtained from 
equations 3, 13 as fol lows 



1 L 



l L 


< p’ ” 

* p > lo 

= 0 , 

( p’ ’ 

* P ’ » lo 

= 0 


■ L 



.L 


( Q’ ' ' 

<50 1 lo 

= 0, 

( Q’ ’ 

cs 

JO 

o“ 

= 0 


{ R’ <5R ) | Q = 0, 


3. 23 


3. 3 Finite element formulation. 

As closed form solution of the system of equations 3.22 
is difficult to obtain, one has to apply one of the approximate 
numerical methods such as Gaierkin, Collocation, FDM, FEM etc to 
get the approximate solution. Since the results obtain by FEM 
converge to the exact solution with the refinement of the mesh, FEM 
is utilized to solve the system of equations. 

The beam is descretized into finite number of two noded 
elements < Fig 3. 2. a). There are five degrees of freedom at each 
node. Over a typical element shown in Fig 3.2. b, { X ) is 


approximated as 
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{ X* } T 


where , 


C S ] = 


0 

0 


t S ) t * ns ) 


N 1 N 2 0 


0 

0 0 


0 


3 


N 1 N 2 0 


M 0 


0 

0 


. . . 3.24 


0 

4 0 

O M , 


and 


{ X n ® > 


P P 0 0 

1 1 *1 


R 1 P 2 P 2 


Q, 


R„ 


... 3.25 

Her© N.( i=l,to 4) and M. (i =1, 2)are the shape functions of the 
local co-ordinate For th© e^ element the relation between £ and 
the global co-ordinate x is given by < Fig 3.2.b ) 


( x - x ) 
? = O 


3. 26 


where. 


{ x , - x ) 

e + 1 e 


( x + x ) 
e+ 1 e 


3. 27 


The expressions for N.’s and M.'s which satisfy the convergence 

r I 1 
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criteria are 

*1 * r 1 2 - 3( * ( 3| 

n 2 » 

N 3 = ^-(2*3f-f 3 ) 

n a * it < - 1 - « * ? 2 * <r 3 I 

"i = 5- ' 1 - < ’ 

M 2 “ 5- < 1 « f 1 


* * * 3.28 


Substituting approximation 3.24 in equation 3.. 21, we get the 
following residue for a typical element 

cx » • • * » » ne 

{*?} = {[ A 3£ S 3 + C B 3 [ S 3 + C C 3C S 3 ) { X ) 

2 no 

- X. £ D ][ S ]{ X J 

. . . 3.29 

Where, £ ) is a non-zero residue vector. To minimise the 

residue, we use the galerkin method. Then we set 


S e+1 C S ] T £ B dx = 0 ... 3. 30 

X 

e 

After performing the integration by parts we get the following 
equat ions 


C K® 

3 { 

x ne ) 


- Ol) 

[ M® 3£ 

X ne 3 = { 

R e ) 

where. 









0) 

II 

y 

fO 

i s 

the 

eigen value. 



C 

K® 3 

is 

the 

e 1 ement 

stiffness 

matr i x , 


C 

M 9 3 

is 

the 

e 1 ement 

mass matrix. 
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sufici i R ) is the element force vector. 

The expressions for these are as follows : 


[ K 


e 


X 



T 1 

{ [ S 3 


+ [ 


C A 3 C S 3 - C S 

3CC3CS])dx 


3 ' C B 3 t S 3 ’ 


C M® 3 = f e+1 { c S T 3 C D 3 [ S 3 3 dx 

x 

0 


, . . 3.32 

. . . 3.33 


C R S 3 


i C S 


T 3 C A 3 C S 3 ’- E S T 3 C A 3 C 
C S T 3 C B 3 C s’ 3 3 { x n ® 3 I 


« « t 

S 3 


x 


0+ 1 


X 

e 


. . . 3 * 34 


To evaluate the element stiffness and mass matrices, we change 
the variable from the global co-ordinate x to the local 
co-ordinate? using equation 3.26. We see that 

... 3. 35 


to -1 and 1 

f<„) dx . L / 1 f( x o * L.f » « _ _ 3 36 

X ~ 1 

e 

Using relations given in 3.35 and 3.36, the integrals 3.32 to 


dx 

d 

dx 


L d£ 
e N 


d 

*K 


and 


dx 


diC 


and the limits of integration change from x^ and k q4 1 
r espec t i v e l y . Thus a typical integral transforms to 
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3.34 are first transformed and then integrated numerically 
using Gauss - Legendre integration scheme to find the stiffness 
and the mass matrices. The element matrices thus obtained, are 
assembled together to form the global stiffness and global mass 
ma t r ices . 

For the beam under consideration, the end x^ = 0 is fixed while the 
end x ^ = L is free. Then the geometric boundary conditions are 
P = 0, P* = 0, Q = 0, Q.’ = 0, R = 0 at x 1 = 0 

.... 3. 37 . a 

Then the < natural ) conditions at x^ = L can be found from 
equations 3.23. They are 

P = 0, p' ' = 0, q’ ’ ’ =0, q’ = 0, R = 0 . . . 3. 37 . b 

When there is an hinge support at an intermediate location, the 
additional geometric conditions at this support are 

P = 0, Q. = 0, R= 0 ... 3. 37. c 

0 

When the element force vectors { R ) are assembled, the 
contributions to an interior node from the neighboring elements 
cancel each other. The contributions to the node at x = L become 
zero because of the boundary conditions 3.37.b. At the node ~0 

the global force vector is unknown. However, when the geometric 
conditions are applied, the number of unknowns get reduced and then 
these equations with unknown right hand side can be discarded. 

After application of the boundary conditions we get a eigen 
value problem of the form 
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C K g 1 { X g } = “ t M g ] { X p } ,,.3.38 

where [ anc * C are modified global stiffness and 

assemb 1 ed ma a> o ins, t r i x r 0 s p?c t i v 8 1 y and { X } is fcH© rna t r i x 

& 

O 

containing the global degrees of freedom. The natural frequencies of 
the system are obtained as square - root of the eigen values of the 
prob l eni 3, 38. 


3. 4 Solution 

The beam under consideration { i.e the VGVAWT blade with 

uniform cross - section ) is shown in Fig 3.3 along with the 

necessary details. We observe that the stiffness matrix C K 1 and 

g 

the mass matrix [Ml are unsymmetric matrices. To solve the 

g 

equation 3.38 we use the NAG subroutine F02BJF. As the continuous 
system has been approximated by a set of finite elements, the 

infinite degrees of freedom corresponding to the actual system 
have been reduced to a finite number. Hence the approximate 

model is more stiff than the actual system. Obviously with the 
increase in the number of elements, the number of degrees of 
freedom increases, the overal l stiffness of the approximate 

model decreases and consequently the eigen values converge to 
their true values. 

During the computation of the results, it was experienced 
that as the number of elements was increased , the first eigen 
value showed a slight oscillation about the mean value, wher eas 
the higher eigen values converged. It was conjectured that this 
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could be possibly due to 


1)1 > > l 

z i z i y i y i 


leading to the 1 1 l -cond i t i on i ng of th 


mat r ices , 

2) The number of iterations in the NAG subroutine are fixed 
but unknown and the user cannot exercise any control over it. The 
eigen values are the roots of the polynomial given by 


C K 3 ~ w C (1 3 I 

1 g g 5 


= 0 


, 3,39 


Th© higher roots of the polynomial equation 3,39 can be obtained by 
neglecting the lower degree terms. But for finding lower eigen 
values particularly the lowest one, one has to consider all the 
terms. Consequently, higher roots can be computed with fewer number 
of iterations, while the first few roots are obtained after a 
considerable number of iterations. Therefore it is possible that 
the first eigen value does not converge within the number of 
iterations fixed by the subroutine i tse l f . F i na l l y it was decided to 
choose 33 elements ( Fig 3, 2, a ), 

3.5 Results arid discussion 

To check the validity of the finite element model, the 
solution pattern obtained from the same is compared with the one 
given in reference C d ] for a beam hinged at both the ends. In our 
case, the cross — section of the beam has an airfoil shape as shown 
in Fig 3*3 while in reference C 6 ] it has the shape of a channel* 

Thus both have one axis of symmetry. Further the geometric and the 

different in these two cases. Figure 


material properties are also 
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3.4 shows the validation of the first three natural 
with Cl, These frequencies are 

1) first flap-wise bending ( bending about 
) First edge-wise bending { bending about 
3 > first torsion 


f requencies 

- axis ) 
axis 1 


It is seen from Fig 3,4 that the edge-wise bending frequency 
increases with the increase in the rotational speed. The first 
flap-wise bending frequency decreases initially with the inorease 
in the speed of rotation until it becomes zero at a particular 
speed. Beyond this speed the system does not vibrate in the first 
flap-wise bending mode. This is indicated by a negative eigen value 
(i.e an imaginary frequency ) beyond a particular speed. This is 
the first critical speed of the system. This bending mode 
of vibration again re-appears after a certain speed of 
rotation ( as indicated by a positive eigen value. ) This is one 
of the higher critical speeds of the system. The frequency 
corresponding to the first torsion mode decreases with the increase 
in the rotational speed and it becomes zero at a particular speed. 

The results of reference C 8 1 are reproduced in Fig 3,5, It 
can be seen that the trends of variation of the frequency with O 
are ident ica l . 

The blade of VGVAWT is a beam fixed at one end and hinged at 
an intermediate point. Before getting the results for the same, it 
was decided to study the trend of variation of frequency with O for 
a cantilever beam. These results are shown in Fig 3.6. It is 

observed that the trends are exactly the same as those for a beam 

the results for the blade of 


hinged at both the ends. Finally 
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VGVAWT ( uniform cross - section ) are shown in Fig 3.7. It is seen 

that the frequencies show the same trend with rotational speed. It 

is observed that, in this case the frequency of a particular mode 
lies in between the corresponding frequencies of a cantilever 
beam and the beam hinged at both the ends, 

A parametric study was conducted to examine the effects of 
wall thickness, shear centre position, the hinge location and taper 
on the natural frequencies. It is observed that the frequencies do 
not change with the wall thickness. This can be explained as 
foil ows : 

Note that al l the geometric properties of the cross - section { 

A, l , I , l and J ) are linear functions of the wall 

y i y i z i 2 i p 

2 2 2 2 

thickness. As a result, the ratios a , b , c and d ( defined by 

equation 3.19 ) are independent of the wall thickness. Further, the 

shear centre location does not change with the wall thickness as it 

depends only on the shape of the cross - section and not on its 

size. Equation 3.20 shows that the dependence of the natural 

frequencies on the geometry of the cross - section is only through 
7 2 2 2 

the ratios a , fa , c , d ( henceforth called the shape factors ) 
and the shear centre location y^ . Since these quantities do not 

depend on the wall thickness, the natural frequencies do not depend 
on the size but only on the shape, provided of course the thickness 
is smal l compared with the other cross - sectional dimensions. 

If an airfoil different from NACA 0015 is used for the cross - 
section, the natural frequencies will change. Calculation of the 
shape factors a 2 , b 2 , c 2 , d 2 and the shear centre location for 

various airfoils is a very tedious task. Therefore, the pa rametr ic 
study with the shape as the parameter is not carried out. Instead, 



only th© effect of unsymmetry of the 


cross 


6 4 


- section on the natural 

frequencies is studied by varying y . This variation is shown in 

k 

Fig 3.8. It is seen that only the torsional frequency changes with 
the shear centre position. It decreases as the distance of the 
shear centre from the centroid increases. For effective functioning 
of the reef control mechanism, the location of the intermediate 
support ( hinge ) should be close to the base v The effect of 
changing the hinge location on the natural frequencies is shown in 
Fig 3*9, It shows that the frequencies increase with the increase 
in the support distance from the base. Since the hinge support 
should be at the mass centre of the beam, a decrease in the support 
distance can be achieved by bringing the mass centre closer to the 
base. This is possible only when there is a taper. One such 
arrangement ( stepped internal taper ) is discussed in chapter 2 
(Fig 2.2). The effect of this taper on the natural frequencies is 
is illustrated in figures 3. 10. a 3.10.b and 3.10.C. It is observed 
that the flap - wise bending and torsional frequencies decrease 
with the taper but the edge * wise bending frequency increases. The 
increase in the edge - wise frequency is due to the^in z^z^ 
over th© tapered portions BC and CD IFig 2.2). 


3.6 Interference diagram. 

The normal and tangential aerodynamic forces acting on the 
blade are periodic in nature. Here we concentrate on the normal 
component since it is the more predominant of the the two 
forces, As seen from fig 2.3 the normal force is an odd function 
of the azimuthal angle and has a period of 2- radians and can be 



represented as b. f our isr 


sine series as 


F (<p) 
n 


oc 




b sin( n# ) 
n 


b s i a ( nft t ) 

n 


where. 


b 

a 



F (<£) sin (n$) d<p 
n 


The coefficients t> n can be found by using a suitabl 
integration scheme. Table 3.1 shows the first 


J, 39 

. . 3.40 

e numerical 
twenty five 
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be small. On the other hand, the margin between the lowest speed of 
interference ( 450 r pm ) and the operational speed ( 85 r pm ) is so 
large, that this small decrease would not cause any resonanae. 

When the ratio of the spar thickness to skin thickness is 
increased, the shear centre moves away from the centroid causing a 
decrease in the first torsional frequency ( Fig 3.8 ). However, 
that would not change the lowest speed of interference (Fig 3.11). 
Thus, again there would be no resonance. 




Cross-section shape — NACA 0015 

c = 0.45 Id, = 0.2 , t x = 6 mrn 

1 = 1.4 m 1 = 3.25m 

Fig 3.3: VGVAWT blade with Uniform Cross-section 
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N-8- The encircled numbers represent 
the element numbers. 

Fig 3. 2. a : Descretization of the Beam into a Finite 
Number of Two-noded Elements. 


x 



2 ,th element 

X -Global co-ordinate, £ -Local co-ordinate 
a -Half the length of the element 
Numbers inscribed in the triangle represent 
local node numbers. 

FIG. 3.2-j?' A TYPICAL TWO-NODED ELEMENT 
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. 10. b : Effect of Taper on the First Edgewise Frequency 







.11 : Interference Diagram for VGVAWT blade. 
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CHAPTER 4 

ANALYSIS OF COMPOSITE BOX SHAPED SPAR 

4* 1 Introduction: 

The wind turbines of high capacity generally rotate at high 

speeds. The blades of such turbines severe loading (both 

aerodynamic and centrifugal forces). In such cases the use of fibre 
wound D-shaped spar may not be sufficient to sustain the 

operating loads. Hence the spar is made of box shaped structure 
consisting of composite laminates stacked over each other. 

Depending on the stacking sequence (layup), there are 
three different cases as shown in Fig 4.1. These are 

1) Symmetric ply layup 

2) Antisymmetric ply layup 

3) cross ply layup 

The present chapter deals with the static analysis of such 

composite box shaped spar, subjected to various types of load such 

as axial force, torque and bending load. In this case there is a 
coupling between different modes of deformation. This coupling 
depends on the structure of each laminate and on the ply stacking 
sequence. The variational finite element formulation is used to 
study the deformation of the spar. 

4* 2 Stiffness coefficients* 

We assume that the spar carries a major portion of the 
load and hence the entire blade cross-section can be represented by 
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the spar for the purpose of static analysis. The spar is modelled 
as a rectangular box beam whose walls consist of laminated 
orthotropic composite plies. 



4«c!:*£:» Constitutive relations: 

The box - beam geometry and co-ordinates are shown in 
Fig 4*2. The origin of the cross-section axes y) - ( is at the 
centroid of the beam* The deformation of the beam is described in 
terms of extension, bending, torsion and transverse shear, The 
effect of constrained warping is neglected. These assumptions 
yield displacements of the form 



In the above equations U, V, and W represent the total axial 
deformation (in the x direction), the total horizontal deformation 
(in the y direction) and the total vertical deformation (in the z 
direction) respectively. u, v and w represent the axial 
deformation due to extension, horizontal deformation due to bending 
and vertical deformation due to bending respectively. <p the twist 
while y represents the transverse shear strain in the cross 
section. 

The rotation of the cross-section about the z axis (ft ^ ) 
and that about the y axis {ft ) are given by 
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ft = v ’ { x ) 
2 

ft = w’ ( X ) 

y 



, , . 4. 2. a 
, , , 4, 2. b 


Substituting 4, 2. a and 4,2,b in equations 4,1. a, 4.1.b 
and 4.1.c, we get the equations for displacements as 


U ( X ) “77 

<X> “ K ft < x ) , 

z y 

. . .4. 3, a 

V { X ) “ £ 

<£ { x ) , 

... 4, 3. b 

w ( x ) + 77 

<^ ( X > . 

... 4, 3. c 


Strains are determined by differentiating these 
displacement terms. Since the walls of the box - beam are assumed 
to be relatively thin, only the axial and in-plane shear strains 
are considered. The resulting strains are 


e = <?U = u* ~ Tj ft - £ ft , ...4. 4. a 

KK #F“ z y 

SU 0W 

£ x% + 

= w’ - + 7) <p’ 

= + r > *' ’ . 4. 4. b 

m + gv 

£ xyj ~ tf/j <?x 

= v ’ - n z - z <t>’ 

- y° - £ ... 4. 4. c 

XT/ 

For a ply in the vertical beam wall, the elastic 


constitutive relations are 
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-- - 1 


— — 

a 

XX 


Q li 

Q 12 

CO 

a? 


e 

XX 

KK 

“ 

^2 

Q 22 

q 26 



i 

<5 

x 

S-* 

1 


^16 

®26 

CO 

CO 

C3T 




« * . 4, 5, a 

For a ply in the horizontal wall the relations are 


— 


1 



<y 

XX 


Q n 

Q 12 

Q 16 


€ 

XX 

cr 

m 

= 

Q 12 

Q 22 

Q '26 


& 

_ XT/ _ 


ft 16 

Q 26 

q 66 


£ 

_ XT/ _ 


• • • 4 * 5 • h 

The terms etc. are the reduced stiffness coefficients of 

a. single ply which are related to the material properties by the 

relation C 9 i 


where 


C T 1 


-{ Q d ) = t T 3 i Q ), 

R 


, . . 4.6 


4 

m 

4 

n 

4 

4 

m 

2 2 
i n 

2 2 
m n 

2 2 
m n 

2 2 
m n 

3 

m n 

3 

-mn 

3 

mn 

-mn 


0 2 2 

2m n 

, 2 2 

4m n 

2 2 

2m n 

2 2 

4m n 


2 

0 2 2 
-2m n 

# 2 2. 

( m - n ) 


4 4 

m + n a 

( m' 3 ~n mn- ) 


2 2 

-4m n , 

mTt 

2(mn 3 - .ranf ) 

' 3 , 


1*3 3 31 "' "J>. 

(jmh - mn ) 2(ra-'Ti - ) 


4.7 
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< > 


{ Q } 


| ^11 11 22 Q 66 ft 12 Q 16 q 26 | 

_l 


, . , 4, 8 


| ®li °'22 Q 12 Q 66 


« « * A * 9 


In expression 4,7, 

m = cos B 9 n = sin B ,,,4,10 

where B is the angle made by the fibres with the x - axis 
( Fig 4,3 ), The elements of { Q. ) are given by 


Q. 


11 


22 


( 1 " V 12 V 21 * 

= V 21 E 1 

< 1 - V 12 V 2l’ 


( 1 " V 12 V 21 » 


ft 66 G 12 


4. 11 


where is the Young’s modulus of the ply in the direction of the 

fibre, E 0 is the Young’s modulus in the direction transverse to the 
fibre, G is the in - plane shear modulus and v and v are 

1 X *£* Zi X 

the in - plane poisson’s ratios. The poisson’s ratio are related by 






4, 12 


Adequate treatment of the transverse in - plane stresses 
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a T )Tl anC * anc * stra * ns and is important in composite 

box - beam analysis. Elastic properties within the composite plies 

change with the ply orientation angle and the poisson’s ratio. The 

anisotropic elastic charecter is t ics of composite plies can result 
# 

in highly two dimensional elastic behaviour. There are three 
methods to account for the in-plane elastic behaviour C9], 

Method 1: 

This is closest to the standard beam theory. Based on the 
initial kinematic assumptions of the beam we assume 

s 0 ...4.13 

mi K * 

and c/ and cr„ H have no effect on the net cross-section force 

nn K < 

and moments. This assumption is very similar in nature to the plane 
strain assumption. This assumption does not take into account the 
in -plane behaviour of the ply very accurately. 


Method 2: 


Here we assume 
= 0 


a 




nn KK 

and £ and are removed from the constitutive 

rrn ? Z 

by substitution. This assumption is similar in nature to 
stress assumption. This is the method used in the present 


4, 14 

relation 
the plane 
anal ys is . 


Method 3; 

This method treats the transverse in-plane behaviour in a 
more refined manner. Condition of in— plane stress and in— plane 
strain are imposed such that there are no in - plane forces or 
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moments , 

Using method 2 and putting = 0 in equation 4,5. a, we get 
the modified constitutive relation for a ply in the vertical wall 


O 1 

XX 



Q ie 


£ 

XX 

1 

X 

1 


Q 16 

^'66 


1 

C!> 

X 

1 


Similarly substituting 0 in equation (4.2.e) we get the 

modified constitutive relations for a ply in the horizontal wall as 


<y 

x x 


0 1 

XT) 


where. 





Q 16 


£ 

X X 

Q '16 «66 


ft 

X 


, . , 4, 15, b 


Q • 
*11 


'll 


( a u» 


Q 


22 


Q ' 
*16 


Q - { Q Q ) 

*16 12 *26 


Q 66 = Q 66 


( Q 26 > 


22 


22 

2 


... 4. 16 


Substituting the relations for strain 4*4. a, 4* 4. b and 

4.4.C in equations 4. 15. a and 4. 15. b, we get the following 

equation for a ply in the vertical wall 
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{ cr } = C V 1 { £ } 
v 


, , , 4 , 1 7 . a 


Similarly for a ply in the horizontal direction we get 


{ cr u } = CH3{£?> 

H 


. . , 4. 17, b 


Here, 


{ & ) 
v 


4. 18 


{ O' } 
n 


. , . 4, 19 


{ £) 


1 u ’ 

o 

o 


ft v 

1 

1 

xr/ 

x< 


y 

2 I 


C VI 


[ H 3 = 


4. 20 


°li 0 °16 n tt 16 


? Q l ’ 1 - r, Vi 


a L 6 0 «66 " °66 - « 0 16 - 71 V 

. . . 4.21 


Q 16 ° " * ft 16 - * Vl ~ 77 Vl 


Q 66 ° ~ ? Q 86 * Q 16 77 Q 16 


. .. 4.22 
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The net forces and moments acting over the cross-section are 
defined to the stresses in the beam walls by equilibrium as 
foil ows , 


F 

. X 


/ / ct dA 

XX 

F 

y 

= 

S S a dA 

F 

z 

= 

/ / o* „ d A 

X? 

fl 

X 

as 

/ / { r) & xl , - K 

M 

y 

= 

- f f <y £ dA 

X X 

M 

z 

= 

- f f <y n dA 

XX 


. . . 4, 23 


where, F , F and F represent the net force in the x, y and z 
x y z 

directions respectively, while M , M and M are the moments 

x y z 

about the x-axis , y-axis and z- ax is respectively* Substituting 
the relations 4. 17. a. and 4.17.b into 4.23, we get 


{ F ) = C D ] { e ) 

where, 


1 F 

F 

F 

M 

M 

M 1 

L * 

y 

z 

X 

y 

2 _l 


. . . 4.24 


... 4 * 


25 


and 


D 1 

i s 

the symmetric 

stiffness matrix whose elements are given 

D n 

= 

/ 

H 

X 

Vi 

dA 

+ S J Q ’ dA 

V 

D 12 

= 

/ 

H 

X 

*16 

dA 


D 13 

= 

f 

V 

f 

*16 

dA 



D 


14 


/ f 
V 



dA 


/ S 
H 


16 


dA 
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15 


/ / * Vi dA 


/ / ? Q,i dA 
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16 


f f 0 Q 1 ' 1 dA 


Sir] dA 


>22 * , / ' Q 66 dA 


°23 = ° 


D 24 ' \ ' ( °66 dA 
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°25 

H 

/ / 

K 

*16 

dA 

°26 
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/ 

H 

f T) 

*16 

dA 

D 

= S 

/ G ’ 

dA 


33 
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66 



°34 

= s s 

V 

17 

*66 

dA 

°35 

V 

f s 

17 


dA 

°36 

V 

S f 

? 

*16 

dA 

°44 

V 

f f 

n 

2 1 

66 

°45 

H 

/ / 

K : 

~ Q ’ 

W 16 

< 

°46 

H 

f s 

77 K 

*16 

dA 


l 66 


s s TJ * dA 
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>55 " ' / * Q U dA 

H 


X X K dA 


>56 * „ ' ' » < Q,’ dA 


f s n £ Q ± ^ dA 


XX Y! Q \ t d A 


xx Tf Q’ il d A 


. . . 4.26 

In the above equations, the suffix H and V represent the 

integration over the horizontal and vertical wall respectively. The 
integrals in 4.26 can be numerically integrated to find the 
stiffness coefficients. For the present case of box-beam, it is 
poss ib i e a 1 develope analytical relations for the stiffness 
coef f i c ients. 

4.3 Finite eleiraent formulation 

The beam is divided into a number of elements. The strain 
energy of a typical e th element can be expressed as 


e+ 1 t 

/ is } CD) is 


} dx 


4.27 


where [ D 3 and is .) are given by equations 4. band 4. 3. h 

e 

respec t i ve l y , Us i ng relations 4.1.d and 4. l.e the matrix { £ ) 

can be expressed as 

{ ^ } = { P ) + i Q. ) . . , 4, 28 


where. 


{ P > 


ti V 


% (¥ ft' 

v z 


{ Q. > 


1 0 

- ft 

- 

0 

0 

0 I 

L 

z 

y 



_I 


... 4. 29 
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Thus, we have 6 degrees of freedom, namely u, v, w, <p, ft and ft . 

y 2 

These can be expressed in the vector form as 


1 u 

V 

w 



/? 1 

1 




y 

£_l 


,.,4.30 


We notice that in equation 4,26 the highest order of derivative 

of any degree of freedom is two. Thus, for convergence, 

\ 

completeness in 1) u, u’ 2) v, v*, 3) w, w* , 4) <£, <p’ , 5) 

and 6) ft , ft* is required and the compatibility of u, v, w, <£, /? 

z z y 

and f3^is Hence two noded elements with linear shape 

functions can be used. The vector { X } at any point can be 


expressed as 




.,,4.33 
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is the vector of elemental degrees of freedom. Here N and N are the 

X a£ 

shape functions and can be expressed in terms of the local 

co-ordinate ou For the e th element, the relationship between oi and 

the global co-ordinate x is given by 

x - x 


a = 


, ,, 4,34 


where, 


L = 


e+ 1 


, 4, 35 


is the half length of the element and 



, , , 4,36 


is the global co-ordinate of the mid-point of the element e. 

The expressions for N and which satisfy the convergence 

criteria are 


N 


1 


1_ 

Z 


1 - a ) 


N 


2 



1 + a ) 


, 4,37 


From expressions 4,28 and 4,29 the vector { e ) can be expressed 
as 

T 

{ £ ) = C B 3 { X n ® ) ... 4.38 

where 


C B ] = C B ] +CB3^ 

2 . ^ 


. . . 4.39 
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0 
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0 0 
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0 -N, 


■n 2 0 


0 o 

o o 

O 0 


. 4.41 


where the subscript , ot represents differentiation w.r.t ct. 
Substituting in 4.27, we get the expression for the strain energy 
as 


4— { X ne ) [ K® 3 { X ns ) 


... 4. 42 
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t K® ] = / C B ) T t D 3 C B ] dx ...4.43 

X 

e 

is the elementary stiffness matrix. To perform the above 
integration numerically, the variable x needs to be transformed to 
the local variable a using the relation 4.34* 

e 

When the elemental stiffness matrices l K 1 are assembled, 
we get the global stiffness matrix C K ]. The global finite element 
equations become 

C K ] { X > = i F } , , , 4, 44 

g g g 


where, { X ) is the vector of global degrees of freedom and { F ) 
is the global force vector. Here it is assumed that the external 
forces consist of only the concentrated forces and/or moments. A 
node is placed at the point of application of every concentrated 


force or moment. Thus { } is found by putting the values of the 

external forces and moments in the approximate rows. 


4.2 Results and Discussion 

The analysis was carried out for a cantilever beam. The 
response of the structure was studied for three layups of the plies 
namely cross-ply layup, antisymmetric layup and symmetric layup. 
The following three types of loading were considered: 

a) Unit tip axial force 

b) Unit tip bending force 

c) Unit tip torque 

The details of the beam are shown in table 4.1 ( Ref. C41 ) 
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a * ^ ross - Ply layup : 

Table 4.2 shows the tip deflections for a cross - ply 
layup. As expected, the various modes of deformation are uncoupled. 
To study the effects of constrained warping, the analysis was 
carried for a beam of having a L/d ratio 29. In this case the 
values of c and d chosen for this case are 52.324 mm and 26.035 mm 



Length of the beam { L ) 


762.00 mm 


Outer depth of the box - beam 


13.64 mm 


cross - section { d ) 


Outer width of the box - beam 


24. 20 mm 


cross - section ( c ) 


Ply thickness ( t ) 


0. 127 mm 


Elastic modulus in the direction 


142.045 Gpa 


of the fibre ( E ) 


Elastic modulus in the direction 


9.796 Gpa 


perpendicular to fibre ( E^ ) 


In - plane shear modulus < G _ ) 


Poisson’s ratio in the direction 


of the fibres ( l 


6.001 Gpa 
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Tab I e 4.2: Response of a cross - ply composite box - beam 


Load ing 

u ( mm ) 

v ( mm ) 

$ ( mm ) 

a) Unit tip axial 

0. 000807 

0.0 

0 . 0 

force 




b) unit tip bending 

0 * 0 

4.072 

0.0 

force 




c) Unit tip torque 

0 . 0 

0 . 0 

0.00334 


Tab l e 4.3 Tip twist under unit tip torque for 

cross - ply layup composite box - beam 


Case 

twist ( radian 1 

without warping 

0.0003749 

with warping 

( ref 1 4 3 ) 

0.000425 


b ) symmetric layup 

Figure 4*4 shows the variation of tip twist with the ply 

orientation angle ( B ) under unit tip torque. It is observed that 

the tip twist reaches its minimum value at S = 30°, The same graph 

also shows the result from reference C 4 ] which indicates the 

effect of constrained warping. When this effect is neglected, the 

tip twist is underestimated up to 6 = 36°and overestimated beyond 

that up to B = 85°, The maximum error in the estimation of the 

twist is about 45 %* When the effect of constrained warping is 

o 

included, the location of minimum twist shifts from & = 30 to 
g? = 45°and the graph becomes symmetric about the 1 ine & - 45 , 
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Figure 4,5 shows the variation of tip twist with the ply 
orientation angle ( B ) under unit tip bending load. It is seen 
that the effect of bending - twist coupling is more predominant 
for the range 15° < 9 < 60°, Further, the coupling effect is quite 
strong as the twist due to unit bending load is almost of the same 
order as that due to unit torque ( Fig 4,4 >, The result of 
reference C 4 ] which includes the effect of constrained warping is 
also shown in Fig 4,5, It is observed that, when the effect of 
constrained warping is neglected, the tip twist is underestimated 
by as much as 70 %, 

Figure 4,6 shows the variation of tip bending deflection 
with the ply orientation angle < B ) under unit tip bending load. 
As expected, the tip deflection increases with S* This happens 
because the effective Young’s modulus along the beam length 
decreases with The same trend is observed in Fig 4,7 which shows 
the variation of the tip extension with B under unit tip axial 
load. In both the graphs, the maximum rate of increase occurs at 
around B = 45°. 

c ) Antisymmetric layup 

Figure 4,8 shows the variation of the tip twist with the 
ply orientation angle ( B ) under unit tip torque. For this case 
also, the minimum value of the tip twist occurs at 9 = 30°, The 
same graph also shows the results from the reference C 4 1 which 
includes the effect of constrained warping. When this effect is 
neglected, the tip twist is underestimated for all values of 6. The 
maximum estimation error is about 25 %. 
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The effect of extension - twist coupling is shown in Fig 

4*9. a in the form of variation of tip twist with & under unit tip 

axial load* The coupling effect is more predominant for 
15 < & < 45°. The result of reference C 4 ] which includes the 
effect of constrained warping is shown in Fig 4* 9* b* It is seen 
that, when this effect is excluded, the tip twist is underestimated 
by as much as 99 %• However, even after including this effect, the 
order of tip twist increases to 0,001 radians which is about 25 % 

of the tip twist due to unit tip torque* Therefore the coupling 
effect is rather weak for the antisymmetric layup. 

Figure 4*10 shows the effect of varying the ply angle ( # ) 
on the tip bending deflection under unit tip bending force* It is 
observed that the bending deflection increases gradually with & 
upto <9 = 60° after which it remains more or less constant* The same 

trend is observed in Fig 4,11 which shows the variation of the tip 

extension under unit tip axial load. 

d ) Compa.r ison 

Comparison of Figures 4*4 and 4.8 show that the tip twist 

under unit tip torque is more or less same for both the symmetric 

and antisymmetric layup* For most values of &, the tip twist for 

symmetric / antisymmetric layup is less than that for the cross-ply 

layup ( Table 4,2 )* The tip bending deflection of the antisymetric 

layup is always smaller than that of the symmetric layup except for 

Q = Q ( Figures 4,6 and 4, 10 1* For & = 90 , it is less than 

one-third of the value of the symmetric layup. Further, the maximum 

tip bending deflection of the antisymmetric layup is less than the 

tip bending deflection of the cross - ply layup < refer to table 
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4.2 K Similar trend is observed for tip extension ( see Figures 
4.7, 4. 11 and table 4,2 ). Thus for the same material and geometric 
properties, the antisymmetric layup is more stiff. 

As stated earlier, the additional twist produced by the 
coupling effect is much less in antisymmetric layup < Pig 4.9 ) 
than in the symmetric layup { Fig 4.5 ), Thus, the coupling effect 
is weak in antisymmetric layup. It is totally absent in the 
cross-ply layup. 

When the effect of constrained warping is neglected, the 
maximum error in the estimation of the tip twist under unit tip 
torque is more for symmetric layup than for antisymmetric layup. 
The total twist consists of the twist due to the torque and the 
additional twist produced by coupling effect. As stated earlier, in 
the symmetric layup, this additional twist is almost of the same 
order as the twist due to torque. On the other hand, it is vary 
small in the antisymmetric layup. Further, when the effect of 


constrained warping is 

i nc l uded. 

the additional 

twist 

in 

the 

symme t r ic layup, becomes 

much larger 

than the the 

twist 

due 

to 


torque ( Figures 4.4, and 4.5 ). But in the case of antisymmetric 
layup, the additional twist even after considering the effect of 
constrained warping, never becomes even as large the twist due to 
torque < Figures 4,8 and 4.9,b ). Thus the effect of constrained 
warping is less severe for antisymmetric layup. 

c ) Cone 1 us ions : 

1) The antisymmetric layup is more stiff than the symmetric 
layup. Additionally, the coupling effect is quite weak in the 
antisymmetric layup. Further, in a typical turbine blade, the axial 
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forces are quite small compared to transverse forces. Then the 
coupling does not have any effect. Finally the effect of 
constrained warping is also less in the antisymmetric layup 
compared to symmetric layup. As a result the antisymmetric layup is 
better than the symmetric layup. 


2) Tip twist 

is maximum 

at 

about 

9 

= 30° 

in 

the 

ant i sy mmetr ic l ay up. 

Further, tip 

bend i ng 

def l ec t ion 

and 

tip 

extension are smaller 

upto 9 = 45°. Thus 

a suitable 

range 

for 

ply 


o o 

orientation is from 9 = 15 to 9 = 45 . 

3 ) For the antisymmetric layup, neglecting the effect of 
constrained warping may not led to serious errors either in values 


or in trends in stress / vibration analysis of wind turbine blades 
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Symmetric Layup Beayn 
* Bending-Torsion Coupling 
- Extension-Shear Coupling 



Anti- Symmetric Layup Beam 
. • Extension-Torsion Coupling 
• Bending-Shear Coupling 


Cross-Pty Layup Beam 

* Alternating Piles at 0 and vu deg 

• No Elastic Couplings 


Fig 4.1 : Composite layup designations. 



Fig 4.2 Box-beam configuration and coordinates. 



Fig 4.3 : Illustration of Ply Oreintation Angle. 







Stacking sequence : Top & Bottom (0) 



(uefpea) aiBue istai dx.T. 


Fig 4.4: Tip Twist of Symmetric Layup Under Unit Tip Torque 




0020*0 



(vretp^j) ai6u« q.sTA^. dxj, 




Tip extension u (nun) 



Ply orientation angle (degrees) 

Fig 4.6: Tip Deflection of Symmetric Layup 
Under Unit Tip Bending Load. 



Fig 4.7: Tip Extension of Symmetric Layup 


Under Unit Tip Axial Load. 
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ply orientation angle (degrees) 

Fig 4.9 Tip Twist of Antisymmetric Layup 
Under Unit Tip Axial Load. 



ply orientation angle (degrees) 

Fig 4.9.t> Tip Twist of Antisymmetric Layup 
Under Unit Tip Axial Load. 





Fig 4.10: Tip Deflection of Antisymmetric Layup 
Under Unit Tip Bending Load. 



Fig 4.11: Tip Extension of Antisymmetric Layup 
Under Unit Tip Axial Load. 



CHAPTER 5 


SUMMARY OF RESULTS AND CONCLUSIONS 

For the VGVAWT blade, a standard 2 cell structure consisting 
of a D-shaped thick-walled spar and thin walled skin is chosen, A 
stepped internal taper is introduced to bring down the mass centre 
of the blade. Uni - directional fabric is used as the blade 
material. The stresses at the critical section are determined 
using Engineer’s theory of bending and Bredt - Batho theory of 
torsion of multi - cell sections. Reasonable values are chosen for 
various thickness ratios. Finally the spar thickness in the 
untapered portion is determined using Tsai - Wu failure criterion 
for composite materials. 


The finite element 

method 

is 

used 

to obtain the 

natural 

frequencies and mode 

shapes 

of 

the 

VGVAWT. The 

periodic 


aerodynamic force acting on the blade is resolved into Fourier 
components to determine the most dominant excitation frequencies. 
The Campbell { or interference ) diagram is constructed to check 
for resonance. It is shown that there is no resonance at the 
operational speed and hence the structural design of VGVAWT blade 
is satisfactory from resonance point of view also. Further,. if a 
small change is made either in structural details ti.e various 
thickness ratios > or in the distance of the intermediate support, 
then also the structure Is safe against resonance. Parametric 
study of the. natural frequencies lead to the following conclusions 

1) In a thin walled beam, the natural frequencies do not 
change with thickness. 

2) The natural frequencies under cons iderat ion increase with 



the increase in the distance of the support from the base. 

3) As the shear centre distance from the mass centre 
increases, the torsion frequency decreases. There is very little 
effect on the flap - wise and edge - wise bending frequencies. 

When a wind turbine is to be designed for high capacity, the 
spar has to be modelled as a composite box - beam. The static 
behaviour of this beam is studied for three different stacking 
sequences : cross - ply layup, antisymmetric layup and symmetric 
layup. Shear deformation and appropriate coupling effects are 
taken into account. The study leads to the following conclusions 

1) The antisymetr ic layup is more stiff than the other two 
layups. Compared to the symmetric layup, the effect of constrained 
warping and the coupling effect are small in the antisymmetric 
layup. Thus the antisymmetric layup is better than the symmetric 
layup. 2) For the antisymmetric layup, the ply orientation angle 
should be between 15°and 45°to minimise the deflections. 

3) If the effect of constrained warping is neglected, the 
stress / vibration analysis of an t i sy mme t r ic layup are still 
expected to give reasonable values / trends. 

Suggestions for future work 

1. The effects of shear deformation and constrained warping 
should be included in the stress analysis of the VGVAWT blade. 

2. a) The effects of shear deformation and constrained 
warping should be included in the vibration analysis of the VGVAWT 
b l ade , 

b) The flexibility of the cross - arm and the tower should 
be included in the vibration study. 
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c) The effects of various airfoil shapes ( NACA and NASA 
series ) on the natural frequencies should be studied, 

3) The effects of constrained warping, stacking sequence and 
ply orientation angle on the natural frequencies of rotating 
composite box - beam should be studied. 
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